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Abstract 

In  this  project,  we  used  computational  models  to  analyze  how  the  intrinsic  dynamical  prop¬ 
erties  of  neural  and  mechanical  systems  interact  to  produce  stable,  but  adaptable  locomotion. 
Animal  locomotion  is  a  rhythmic  behavior  that  requires  the  effective  coupling  of  multiple  feed¬ 
back  loops,  including  mechanical  coupling  between  the  animal’s  body  and  the  environment, 
coupling  between  muscular  force  production  and  body  movement,  and  sensory  feedback.  Flo- 
quet  theory,  a  branch  of  nonlinear  dynamics,  includes  ways  to  analyze  how  such  rhythmic 
systems  respond  to  perturbations.  We  developed  several  robust  ways  of  estimating  the  Floquet 
modes  of  a  rhythmic  system,  which  are  canonical  patterns  of  activity  after  a  perturbation.  We 
found  that  when  a  block  of  muscle  is  forced  to  change  length  sinusoidally  and  is  cyclically 
activated,  it  is  strongly  self-stabilizing,  even  with  no  sensory  feedback.  When  two  muscles  act 
antagonistically,  as  they  do  around  most  vertebrate  joints,  then  the  system  is  less  stable  natu¬ 
rally.  However,  with  sensory  feedback,  the  joint  can  be  stabilized  very  easily.  This  research 
may  be  extended  to  analyze  Floquet  modes  based  on  empirical  data,  to  examine  the  stability 
properties  of  real  muscle,  and  to  study  the  stability  of  fish  swimming  and  control  potential  of 
fish  fins. 

Animals  must  move  effectively  through  complex,  unpredictable  environments.  As  they  move, 
neural  circuits  called  central  pattern  generators  (CPGs)  produce  a  basic  pattern  of  muscle  activity 
(Grillner,  2003),  activating  muscles  that  move  the  body,  which  interacts  with  the  external  environ¬ 
ment  (Tytell  et  al.,  2011).  At  the  same  time,  the  environment  produces  forces  back  on  the  body, 
influencing  the  motion  (Jordan,  1996),  and  the  CPG  receives  sensory  inputs  that  modulate  the  loco¬ 
motor  pattern  (Rossignol  et  al.,  2006)  (see  Fig.  1).  Each  of  these  components  has  its  own  intrinsic 
dynamical  properties,  but  all  of  them  must  work  together  to  produce  a  stable  pattern. 

In  this  project,  we  investigated  how  neural  and  mechanical  systems  can  work  together  to  pro¬ 
duce  stable  cyclical  dynamics.  Specifically,  we  examined  how  such  systems  respond  to  pertur¬ 
bations:  small  disturbances  away  from  the  normal  pattern,  called  a  limit  cycle.  For  some  pertur¬ 
bations,  the  system  rapidly  returns  to  its  normal  limit  cycle;  for  other  perturbations,  the  activity 
pattern  changes  for  long  periods  of  time. 
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A  branch  of  dynamical  systems  theory  called  Floquet  analysis  (Floquet,  1883;  Guckenheimer 
and  Holmes,  1983)  allows  one  to  classify  these  types  of  perturbations.  The  state  of  the  system  x(t) 
can  be  defined  in  a  high-dimensional  state  space  (Mp),  where  x  might  include  all  of  the  membrane 
potentials  and  synaptic  conductances  of  neurons  in  the  CPG  model,  or  the  lengths,  velocities,  and 
calcium  concentrations  in  a  muscle  model.  Since  the  system  is  periodic,  x(t)  traces  out  a  closed 
loop  7  in  the  n— dimensional  space.  A  perturbation  means  that  the  system  will  deviate  from  the 
limit  cycle: 

z(t)  =  x(t)  -x*(t),  (1) 

where  x*  (t)  is  on  the  limit  cycle  7. 

For  a  particular  system,  Floquet  analysis  allows  one  to  define  a  set  of  canonical  modes  u k(t) 
that  are  particular  patterns  of  deviations  from  the  limit  cycle.  The  modes  are  the  axes  of  a  new  co¬ 
ordinate  system,  centered  on  the  limit  cycle,  that  rotates  and  stretches  about  the  cycle.  Even  though 
the  modes  may  be  complicated,  when  the  state  of  the  system  is  expressed  in  this  new  coordinate 
system,  the  dynamics  become  very  simple:  any  perturbation  simply  decays  exponentially  back  to 
the  limit  cycle.  The  rate  of  decay  is  quantified  by  the  Floquet  exponent  fj,k,  so  that  the  deviation 
from  the  limit  cycle  can  be  written  as 

z  (t)  =  e^ktu  fc(t).  (2) 

Negative  Floquet  exponents  correspond  to  stable  patterns,  and  the  more  negative  the  exponent,  the 
more  rapidly  the  perturbation  dies  away. 

The  goal  of  this  project  was  to  examine  how  the  intrinsic  dynamics  of  neural  and  mechanical 
systems  interact  in  a  computational  model,  based  on  the  lamprey.  The  primary  results  of  the  work 
were 

•  Several  robust  techniques  were  developed  for  estimating  limit  cycles  and  Floquet  modes.  In 
particular,  we  developed  a  technique  for  finding  Floquet  modes  that  does  not  require  any  inte¬ 
gration.  These  techniques  may  be  used  when  the  state  equation  is  known,  but  also  potentially 
on  empirical  data. 

•  We  found  that  muscle  is  strongly  self-stabilizing  when  activated  cyclically,  possibly  because  of 
the  nonlinearity  in  how  calcium  binds  and  releases  from  muscle  filaments.  Two  muscles  in  an 


CPG  muscle  body  environment 


proprioception 


Figure  1:  Diagram  of  the  nested  feedback  loops  that  make  up  a  neuromechanical  system.  A  central 
pattern  generator  circuit  activates  muscle  that  produces  force  to  bend  the  body,  which  interacts  with 
the  environment.  At  the  same  time,  the  environment  produces  reaction  forces  back  on  the  body 
(“mechanics”),  muscle  force  is  coupled  with  body  motion  due  to  length  and  velocity  dependence,  and 
the  CPG  receives  proprioceptive  sensory  inputs. 
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antagonistic  arrangement,  similar  to  that  of  most  vertebrate  joints,  have  a  less  stable  mode  that 
would  need  to  be  corrected  by  sensory  feedback,  but  may  also  provide  a  simple  way  to  control 
turning. 

•  Preliminary  evidence  was  found  for  an  optimal  feedback  strength  for  maximum  stability  in  the 
neuromechanical  system. 

These  findings  suggest  several  avenues  to  continue  and  extend  this  work: 

•  The  harmonic  balance  technique  may  be  appropriate  for  extracting  limit  cycles  and  Floquet 
modes  from  data.  Further  work  will  be  necessary  to  establish  the  feasibility  and  accuracy  of 
this  approach. 

•  The  analysis  of  the  muscle  model  makes  some  strong  predictions  about  the  stability  of  cycli¬ 
cally  driven  muscle.  The  model  approximates  a  standard  experimental  configuration  called  a 
“work  loop”;  thus,  these  predictions  may  be  possible  to  test  experimentally. 

•  The  ideas  developed  here  may  be  applicable  to  understanding  the  stability  and  control  of  fish 
swimming.  By  using  stimulating  electrodes  to  introduce  phase-locked  perturbations  to  the 
muscles  of  a  swimming  fish,  we  may  be  able  to  rigorously  evaluate  the  phase-dependent  control 
potential  of  fish  fins,  and  how  fish  use  them  to  maintain  stability. 

Details  of  both  the  findings  and  the  future  directions  follow. 


1  Primary  findings 

1.1  Robust  techniques  for  estimating  limit  cycles  and  Floquet  modes 

Much  of  the  project  period  was  spent  developing  numerical  methods  for  robust  estimation  of  limit 
cycles  and  Floquet  modes.  The  numerical  problem  is  made  particularly  difficult  because  we  needed 
to  detect  relatively  tiny  deviations  from  the  relatively  large  limit  cycle. 

1.1.1  Limit  cycles 

Shooting  methods  Initially,  we  proceeded  based  on  the  idea  of  a  Poincare  section:  a  hyperplane 
of  co-dimension  P  —  1  that  intersects  the  limit  cycle  transversely.  The  (linearized)  return  map,  J, 
(also  called  a  Poincare  map)  describes  locally  how  perturbations  from  a  point  x*  (t)  recover  after 
one  full  cycle  of  the  nonlinear  dynamics.  Given  a  system  defined  as 

x  =  /(x),  xef,  (3) 

one  can  integrate  the  dynamics  of  a  point  x0  that  starts  on  the  hyperplane  until  it  returns  to  the 
same  plane  at  another  point  xi. 

As  a  method  for  estimating  the  limit  cycle,  this  forward  integration,  called  a  “shooting”  tech¬ 
nique,  works  well.  One  starts  with  a  candidate  point  close  to  the  limit  cycle  x/P)0  and  integrates  to 
find  the  first  return  xjpl.  Then,  using  a  standard  minimization  technique  such  as  Newton’s  method, 
one  can  minimize  the  difference  |xyp  0  —  x jp.  1 1  • 
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Harmonic  balance  technique  One  can  also  estimate  the  limit  cycle  as  a  Fourier  series.  This 
method  is  attractive  because  it  implicitly  enforces  the  limit  cycle  solution  to  be  periodic.  Full 
details  are  given  in  Appendix  A.  The  approximation  method  is  quite  similar  to  that  used  by  Bo- 
nani  and  Gilli  (1999).  For  system  (3),  choose  the  coefficients  of  a  Fourier  series  that  produces  a 
candidate  cycle  x.  If  x  were  on  the  true  limit  cycle, 

£  -  /(x)  =  0.  (4) 

The  derivative  of  a  Fourier  series  can  be  easily  expressed  in  terms  of  the  coefficients  of  the  original 
series,  which  allows  one  to  simplify  x.  Furthermore,  (4)  holds  at  each  point  along  the  limit  cycle, 
which  means  one  can  discretize  the  problem  into  M  equations,  each  at  a  different  phase.  Then  the 
coefficients  of  the  Fourier  series  can  be  found  by  standard  multivariate  minimization  procedures, 
such  as  Matlab’s  f  solve  function.  See  Appendix  A  for  details.  Interestingly,  this  procedure  could 
also  be  used  to  find  an  unstable  limit  cycle. 

1.1.2  Floquet  modes 

In  principle,  one  can  extend  the  shooting  technique  to  estimate  the  linearized  return  map  J(T). 
J(0)  =  Ip,  the  P-dimensional  identity  matrix.  J(f)  can  be  estimated  by  integrating 

j(i)  =  A(t)J(f),  (5) 

where  A(t)  is  the  derivative  matrix 

A(t)  =  *  (6) 

dxi 

evaluated  on  the  limit  cycle.  Then  J(T)  is  estimated  by  integrating  (5)  from  0  to  T,  where  T  is  the 
period  (Apri  et  al.,  2010;  Guckenheimer  and  Holmes,  1983).  The  eigenvalues  and  eigenvectors  of 
J(T)  correspond  to  the  Floquet  multipliers  My .  (where  My,  =  e//fcT)  and  the  modes  at  time  t. 

However,  this  technique  fails  when  the  exponents  are  large  and  negative,  or  when  they  are  close 
to  zero  but  not  exactly  equal  to  zero  (Fairgrieve  and  Jepson,  1991;  Lust,  2001).  When  exponents 
are  large  and  negative,  perturbations  die  out  extremely  rapidly,  and  the  Floquet  modes  are  entirely 
swamped  by  numerical  error  in  the  integration  over  a  full  period.  Similarly,  it  becomes  very  diffi¬ 
cult  to  distinguish  modes  with  exponents  that  are  exactly  zero  from  those  that  are  close  to  but  not 
equal  to  zero,  for  the  same  reason  that  small  differences  become  swamped  in  numerical  error  over 
the  integration  period. 

A  better  approach  is  to  approximate  the  Floquet  modes  as  Fourier  series.  This  method  is  similar 
to  the  one  developed  by  Traversa  and  Bonani  in  a  series  of  papers  characterizing  the  noise  in  os¬ 
cillatory  electrical  circuits  (Bonani  and  Gilli,  1999;  Traversa  and  Bonani,  201  la,b;  Traversa  et  al., 
2008).  With  some  careful  approximation  and  discretization,  finding  Floquet  modes  becomes  a 
simple  eigenvalue  problem,  requiring  no  integration. 

Linearizing  (3)  about  the  limit  cycle  and  using  (3),  one  obtains 

z  —  A(t)z  =  0.  (7) 

If  ii/,.  is  a  Fourier  approximation  of  a  candidate  Floquet  mode  of  the  system,  and  /4  is  the  approxi¬ 
mation  of  its  Floquet  exponent,  then 

z  (t)  =  eMktu  k(t)  (8) 
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should  be  a  solution  of  (7).  Substituting  into  (7),  we  find 

fikuk(t)  +  ufc(t)  -  A(f)ufc(f)  =  0  (9) 

after  factoring  out  ef‘kt.  Again,  since  we  express  u k(t)  as  a  Fourier  series,  we  can  easily  write  the 
coefficients  of  in  terms  of  the  original  coefficients.  We  discretize  at  M  phase  values  along  the 
limit  cycle.  In  the  end,  the  modes  and  their  exponents  end  up  as  eigenvectors  and  eigenvalues  of 
a  large  matrix  that  is  primarily  constructed  from  the  values  of  the  derivative  matrix  A  evaluated  at 
the  M  phases.  See  Appendix  B  for  details. 

1.2  Intrinsic  stability  of  muscle 

Previous  work  has  suggested  that  muscle  can  be  self- stabilizing  for  static  force  production  (e.g. 
McMahon,  1984),  but  no  one  has  investigated  this  idea  in  the  context  of  cyclical  motions,  which 
are  much  more  common.  In  addition,  most  of  the  previous  work  has  examined  muscle  together  with 
its  sensory  feedback  (e.g.  Prochazka  et  al.,  1997).  We  found  that  muscle,  with  only  its  natural 
calcium  dynamics  and  length,  velocity,  and  tension  relationships,  is  strongly  self-stabilizing.  I 
plan  to  submit  these  results  as  an  abstract  to  present  at  the  Society  for  Integrative  and  Comparative 
Biology  annual  meeting  in  January,  2013,  and  I  anticipate  writing  up  and  publishing  them  in  the 
near  future. 


Brief  description  of  the  muscle  model  The  model  of  muscle  is  taken 
from  Williams  et  al.  (1998),  who  developed  it  for  lamprey  muscle,  but  the 
structure  of  the  model  is  appropriate  for  many  different  vertebrate  mus¬ 
cles.  Briefly,  muscle  consists  of  a  contractile  element  and  a  series  elastic 
element  (Fig.  2).  The  contractile  element  produces  force  Pc  that  depends 
on  its  length  lc,  contraction  velocity  vc,  and  the  amount  Caf  of  calcium 
that  is  bound  to  the  contractile  filaments: 

Pc  =  P0\(lc)Z{vc)Caf,  (10) 

where  P0  is  the  maximum  force.  The  functions  A(/c)  and  £(vc)  define  the 
length-tension  and  velocity-tension  relationships: 


A  (U 

—  i  +  x2(ic  — 

Co) 

(11) 

£(0 

f  f mV ci 

if  vc  <  0, 

(12) 

1 

if  vc  >  0. 

The  calcium  dynamics  can  be  written  with  two  states: 


contractile 


series 

elastic 


Figure  2:  Diagram 
of  the  muscle  model, 
showing  the  contractile 
element  on  the  left  and 
the  series  elastic  ele¬ 
ment  and  damper  on 
the  right. 


Ca  =  (k4Caf  —  k3Ca)(l  —  Caf)  +  Caact,  (13) 

Caf  =  —  {k4Caf  —  k3Ca)(l  —  Caf),  (14) 


where  Caact  depends  on  the  neural  activation  gact,  which  is  scaled  to  be  between  0  and  1: 


Chad  =  h gact(C  -  Ca  -  Caf)  +  k2(  1  -  gact)Ca{C  -  S  -  Ca  -  Caf).  (15) 
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Finally,  the  whole  muscle  element  is  a  spring-mass-damper  system: 

Wlmjc  drn^c  ^ c  ^so)  -^ci  (lb) 

where  M  is  the  mass  of  the  muscle  element  (called  a  sarcomere),  D  is  the  damping,  //  is  the 
stiffness  of  the  series  elastic  element,  L  is  the  total  sarcomere  length,  and  /s0  is  the  resting  length 
of  the  series  elastic  element  (Fig.  2). 

Stability  of  a  cyclically  driven  muscle  segment  Based  on  this  model,  one  can  examine  a  cyclical 
system  in  which  the  muscle  length  is  forced  to  change  sinusoidally  ( L  =  A  sin(27r/f))  and  the 
muscle  is  activated  at  a  particular  phase  <pact  during  the  cycle.  This  is  the  so-called  “work  loop” 
method  used  in  numerous  experimental  studies  (e.g.  Josephson  and  Stokes,  1999). 

For  ease  of  comprehension,  we  use  a  “half-life”  as  a  way  to  compare  the  Floquet  exponents. 
This  is  the  time  that  it  takes  for  half  of  the  perturbation  to  die  out:  ti/2  =  log(0.5)  / ///.,  where  ///;.  is 
the  Floquet  exponent. 

In  the  work  loop  configuration,  the  slowest  Floquet  mode  has  a  half  life  of  10.1%  of  a  cycle, 
regardless  of  the  phase  of  activation.  (Note  that  unlike  autonomous  systems,  this  is  a  driven  system, 
so  there  is  no  mode  with  an  exponent  of  zero.)  All  others  modes  die  away  much  more  quickly, 
in  2%  or  less  of  a  cycle.  Because  of  the  short  time  constant  of  the  mode,  perturbations  have 
a  small  effect  in  general.  Fig.  3a  shows  the  effect  on  force  production  of  a  20%  perturbation 
along  the  first  Floquet  mode.  Note  that  all  effects  are  small  and  die  away  rapidly,  even  though 
the  perturbation  is  very  large  compared  to  physiological  changes,  which  tend  to  be  around  ±5%. 
Similarly,  the  effect  on  total  work  done  is  small,  regardless  of  the  phase  of  activity  or  the  phase 
of  the  perturbation  (Fig.  3b).  This  stability  seems  to  be  related  to  the  calcium  dynamics,  because 
changing  the  calcium  rate  constants  can  increase  the  duration  of  the  first  Floquet  mode.  Increasing 
k4  (the  rate  that  calcium  is  released  from  the  filaments)  or  decreasing  /c3  (the  rate  that  calcium  binds 
to  the  filaments)  tend  to  increase  the  stability  of  the  muscle;  either  change  results  in  less  calcium 
bound  to  the  filaments. 

Stability  of  two  antagonistic  muscles  The  work-loop  configuration  is  simple  and  is  interesting 
because  of  the  large  literature,  but  it  is  not  a  very  realistic  situation.  A  more  realistic  model  would 
be  two  antagonistic  muscles,  pulling  on  a  mass-spring-damper  system,  similar  to  the  muscles  that 
move  a  fish’s  tail  from  side  to  side,  or  to  the  extensor  and  flexor  muscles  that  move  a  leg  forward 
and  back.  Here,  the  two  muscles  are  activated  in  anti-phase  at  a  frequency  fact.  At  the  same 
time,  the  mass-spring-damper  has  a  mechanical  resonant  frequency  fres  and  a  damping  coefficient 
(.  This  configuration  shows  evidence  for  the  necessity  of  sensory  feedback,  but,  at  the  same 
time,  the  possibility  of  using  natural  dynamics  to  steer  an  organism. 

In  this  configuration,  a  new,  slower  mode  appears,  with  a  half  life  of  as  much  as  80%  of  a  cycle 
(Fig.  4a).  The  time  constant  of  the  mode  depends  on  the  mechanical  resonant  frequency  and  the 
damping  coefficient.  The  second  mode  (Fig.  4b)  is  the  same  as  the  first  mode  with  a  single  muscle, 
seen  above,  and  does  not  depend  on  frequency.  The  first  mode  exclusively  affects  the  position 
and  velocity  of  the  mass  (Fig.  4c),  and  the  corresponding  lengths  of  the  muscle  segments.  In  other 
words,  because  two  identical  muscles  pulling  against  each  other  produce  forces  that  are  very  nearly 
matched,  if  the  mass  is  perturbed  off  its  center  position,  then  it  takes  a  long  time  to  return  to  the 
center  position. 
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Figure  3:  Perturbations  have  a  relatively  small  effect  on  a  single,  cyclically  driven  muscle,  (a)  Change 
in  force  due  to  20%  perturbations  along  the  slowest  Floquet  mode.  Black  dashed  line  shows  steady-state 
force  production;  red  lines  show  perturbations,  (b)  Contour  plot  showing  the  change  in  total  work  done 
after  a  20%  perturbation  along  the  slowest  Floquet  mode,  relative  to  the  phase  of  activation  and  the 
phase  of  the  perturbation.  Note  that  all  changes  are  less  than  about  ±10%.  (c)  Half  life  of  the  slowest 
Floquet  mode  as  a  function  of  calcium  parameters.  The  stability  of  muscle  seems  to  depend  primarily 
on  the  calcium  dynamics. 


This  slow  mode  points  to  a  situation  where  neural  and  mechanical  systems  must  work 
together,  each  compensating  for  the  other.  In  this  case,  a  CPG  driving  the  two  antagonistic 
muscles  would  only  need  to  monitor  and  correct  the  mean  position  of  the  mass,  in  order  to  cancel 
out  the  slow  mode  and  produce  extremely  stable  cyclical  motion.  Alternatively,  to  produce  turning, 
the  nervous  system  would  only  need  to  generate  a  brief  perturbation  in  the  position  of  the  mass  to 
cause  a  long  lasting  deviation  to  one  side,  which  would  cause  a  turn  for  a  swimming  organism. 


1.3  Role  of  feedback 


We  are  currently  working  to  simulate  and  understand  the  full  neu¬ 
romechanical  model  shown  in  Fig.  1.  Preliminary  results  suggest 
that  the  system  is  also  quite  stable,  but  that  there  is  an  optimal 
feedback  strength  for  maximal  stability  (Fig.  5).  These  simula¬ 
tions  are  computationally  intensive  to  perform  in  Matlab,  and  may 
require  some  optimization  in  order  to  be  tractable. 


2  Future  directions 


This  project  has  generated  many  useful  ideas  that  should  be  elab¬ 
orated  and  studied  further. 

2.1  Estimation  of  limit  cycles  and  Floquet  modes 
from  data 

A  goal  of  this  project  was  to  develop  numerical  analysis  tech¬ 
niques  that  might  be  appropriate  for  use  on  data,  not  just  on  sim- 
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Figure  5:  Preliminary  results 
showing  the  first  three  non-zero 
Floquet  exponents  for  the  com¬ 
plete  neuromechanical  system, 
shown  in  Fig.  1,  as  the  strength 
of  the  proprioceptive  feedback 
varies.  Note  that  there  seems  to 
be  a  minimum  value  for  the  ex¬ 
ponent  shown  in  blue  at  a  feed¬ 
back  strength  of  2. 


Frequency  difference  f  -  f  (Hz)  Frequency  difference  f  -  f  (Hz) 

act  res  act  res 


(a)  Mode  1 


(b)  Mode  2 


(c)  Effect  of  mode  1 


Figure  4:  With  two  antagonistic  muscles,  a  slower  mode  appears  that  depends  on  the  resonance  and 
damping  of  the  mass  on  which  the  muscles  pull,  (a)  Half-life  of  mode  1  as  a  function  of  the  difference 
of  the  activation  frequency  fact  and  the  mechanical  resonant  frequency  fres ,  and  as  a  function  of  the 
damping  coefficient  (b)  Half-life  of  mode  2  for  the  same  cases  as  in  panel  (a),  (c)  Effect  of  mode  1 
on  the  position  L  and  velocity  V  of  the  mass  for  fact  —  fres  =  0.7  and  £  =  4. 


ulations.  The  harmonic  balance  technique  developed  here  may  be 
adapted  for  estimating  limit  cycles  from  data,  and  may  even  be 
used  for  estimating  Floquet  modes.  I  briefly  describe  how  this 

might  be  possible;  further  work  will  be  necessary  to  establish  the  feasibility  of  this  approach.  As  I 
elaborate  these  ideas  further,  I  plan  to  develop  a  manuscript  for  publication. 

Estimating  limit  cycles  To  estimate  a  limit  cycle  from  data,  one  can  perform  essentially  the 
same  optimization  approach  described  in  Appendix  A.  Given  a  noisy  time  series  data  set  X(tj) 
at  many  times  ti,  one  can  estimate  x  numerically  by  central  differences  as  cZX.  Again,  choose  a 
Fourier  approximation  for  a  candidate  limit  cycle  x.  Interpolate  the  estimated  state  velocity  cZX  on 
to  the  candidate  limit  cycle.  Then,  optimize  the  limit  cycle  to  minimize 

x  -  dX|x=~  .  (17) 

As  before  x  can  be  expressed  in  terms  of  the  original  Fourier  coefficients.  One  would  have  to 
impose  some  smoothness  constraint  on  x;  otherwise,  it  might  wiggle  around  to  every  data  point. 
Also,  one  might  need  to  have  some  outlier  detection  for  velocities.  Some  values  of  c/X  might 
represent  the  true  state  velocity  of  the  system,  but  others  will  be  noise. 

Estimating  Floquet  modes  To  use  the  harmonic  balance  approach  to  estimate  Floquet  modes, 
one  only  needs  to  know  the  limit  cycle  (which  may  be  estimated  using  the  technique  above)  and 
the  value  of  the  derivative  matrix  on  the  limit  cycle.  The  derivative  matrix  does  not  need  to  be 
evaluated  everywhere.  Estimating  a  high-dimensional  derivative  matrix  is  challenging,  but  it  will 
be  significantly  easier  to  estimate  on  a  small  number  of  points,  rather  than  across  the  entire  state 
space.  To  estimate  the  derivative  matrix,  it  is  probably  best  to  introduce  controlled  perturbations, 
rather  than  trying  to  analyze  existing  time  series.  These  ideas  inform  some  of  the  experiments 
proposed  in  §2.3. 
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2.2  Stability  properties  of  muscle 

The  observations  described  above  (§1.2)  make  strong  predictions  about  the  behavior  of  muscle. 
Such  predictions  can  be  tested,  at  least  to  some  extent,  using  in  vitro  muscle  physiological  experi¬ 
ments.  The  work  loop  technique  used  as  a  baseline  model  of  muscle  is  an  experimental  technique 
(e.g.  Josephson,  1993).  Thus,  one  could  perform  physical  experiments  that  are  quite  similar  to  the 
computational  models  developed  here,  in  order  to  test  the  predictions  of  the  model. 

An  experimental  challenge,  however,  is  measuring  and  perturbing  the  system’s  state.  Total 
length  and  velocity  are  straightforward  to  measure  and  manipulate,  but  one  cannot  measure  the 
state  of  the  contractile  element  separately.  However,  one  can  measure  the  elasticity  of  the  se¬ 
ries  elastic  element,  and  to  some  extent  the  passive  damping  properties  of  the  muscle  (McMahon, 
1984),  which  would  enable  an  estimate  of  the  state  of  the  contractile  element.  Calcium  concentra¬ 
tions  are  possible  to  measure  using  calcium  sensitive  fluorescent  dyes  (Vergara  et  al.,  1991),  but 
can  only  be  perturbed  over  long  time  periods  by  changing  total  calcium  concentrations  in  the  saline 
bath  holding  the  muscle  fiber. 

Such  experiments,  while  challenging,  may  help  to  establish  further  the  self-stabilizing  proper¬ 
ties  of  muscle  during  cyclical  motions. 

2.3  Control  potential  of  fish  fins 

Inspired  by  these  computational  results  and  the  experiments  of  Sponberg  et  al.  (2011a,b),  I  am 
developing  experiments  to  examine  the  stability  of  fish  as  they  are  swimming  and  establish  the 
control  potential  of  fish  fins.  In  these  experiments,  we  will  implant  fine  wire  electromyographic 
electrodes  (e.g.  Tytell  and  Lauder,  2002)  in  both  the  axial  red  muscles  that  power  the  side-to-side 
motions  of  the  body  and  tail,  and  in  the  smaller  muscles  that  control  the  motions  of  the  dorsal  and 
anal  fins.  These  electrodes  will  serve  a  dual  purpose:  they  will  both  record  muscle  activity  and  will 
also  be  able  to  stimulate  the  muscle  to  contract.  My  previous  work  has  shown  that  the  dorsal  and 
anal  fins  are  not  passive  keels,  but  actively  generate  substantial  forces  (Tytell,  2006). 

Using  the  techniques  described  above  (§2.1),  we  will  estimate  the  steady  limit  cycle,  including 
both  three-dimensional  kinematic  data  and  muscle  activity.  Then,  we  will  perturb  the  muscle 
activity,  stimulating  the  muscles  at  particular  phases  in  the  swimming  cycle.  This  will  allow  us 
to  measure  the  control  potential  of  the  fins  as  a  function  of  phase,  and  to  estimate  something 
analogous  to  Floquet  modes.  Because  we  cannot  measure  all  of  the  states  of  the  swimming  fish, 
we  will  only  be  able  to  estimate  these  modes  in  a  smaller  subspace. 

Nevertheless,  controlled  perturbations  of  this  sort  will  be  informative  and  perhaps  surprising. 
Sponberg  et  al.  (201  lb)  found  that  a  muscle  in  the  cockroach  leg  produced  a  very  different  effect 
when  it  was  stimulated  during  the  running  cycle,  compared  to  its  effect  during  standing.  Similarly, 
we  might  expect  counterintuitive  results  from  out  experiments.  For  example,  fins  have  muscles 
called  “elevators”  that  elevate  the  fin  rays,  making  the  fin  area  larger.  If  the  fish  is  stationary,  these 
muscles  probably  have  very  little  effect  on  forward  propulsion.  During  swimming,  however,  they 
may  have  a  large  effect,  because  they  change  the  fin  area. 
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Appendices 

A  Harmonic  balance  approach  for  estimating  limit  cycles 

The  approach  described  below  is  quite  similar  to  that  developed  in  papers  by  Traversa  and  Bonani 
(Bonani  and  Gilli,  1999;  Traversa  and  Bonani,  2011a;  Traversa  et  al.,  2008). 

Consider  a  dynamical  system  defined  as 


/(x)  =  0, 


(18) 


where  x  is  a  state  variable  in  P  dimensions.  We  are  aiming  to  write  the  dynamical  system  as  a  set 
of  equations  for  the  Fourier  components  of  the  solution. 

Write  the  limit  cycle  solution  x(f)  in  terms  of  N  harmonic  components  for  each  of  the  i  di¬ 
mensions: 


N 

Xi(t)  =  x{0  +  [xfcik  cos(kut)  +  x{ik  Sin  (hut) 
k=  1 


Based  on  the  coefficients,  define  a  P(2N  +  1)  component  frequency  state  vector 


f  r  /  f  f 
X  —  Lx10  "Gil  xsll 


f  f  f  f  f 

XclN  XslN  X20  Xc21  Xs21 


f  f  T 

xcPN  xsPN J 


(19) 


(20) 


that  consists  of  the  components  for  each  dimension  stacked  on  top  of  each  other. 

Now  let’s  sample  the  system  at  M  samples  in  time,  spread  evenly  through  the  period  T  = 

27 r/cc: 

«=-i_  (21) 

1  u  M  ’ 

where  M  >  2 N  +  1.  Define  a  time-sampled  state  vector,  similar  to  the  frequency  state  vector 

=  [xi(fi)  Xi(t2)  •••  Xi  (tM)  x2(ti)  x2(t2)  •••  xP(tM)]T  (22) 

Then  we  can  write  x*  and  x^  in  terms  of  each  other  by  multiplying  by  a  constant  matrix  T_1: 


the  diagonal. 


where 


=  r^1x/  x/  =  Txk 

(23) 

1  consists  of  an  M  x  2 N  +  1  dimensional  matrix  T0  1  that  is  repeated  P  times  along 

"l 

7cii 

Tall  7c12  7s12 

IclN 

IslN 

r->  — 

1  0  “ 

1 

7c21 

7s21  7c22  7*22 

lc2N 

ls2N 

(24) 

1 

'JcMl 

7sM1  7cM2  7sM2  '  '  ' 

IcMN 

'JsMN 

7 cjk  =  cos(27 rjk/M) 

(25) 

Isjk  =  sin(2vr  jk/M), 

(26) 
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and  finally 


0 

0 


(27) 


V  o 

r-i  =  0  To1 

0  0 

with  P  repetitions  of  Tq  1 . 

We  can  express  x  easily 

N 

Xi(t)  =  E  —x*cikku  sin  [hut)  +  x{ikku  cos  (kut)  (28) 

k=  i 


or 


x-'  =  cefix^, 


(29) 


where  is  a  constant  matrix.  Again,  we  have  an  M  x  2N  +  1  dimensional  matrix  O0 


'o 

0 

0 

0 

0  ••• 

0 

o' 

0 

0 

1 

0 

0  ••• 

0 

0 

0 

-1 

0 

0 

0  ••• 

0 

0 

0 

0 

0 

0 

2  ••• 

0 

0 

0 

0 

0 

-2 

0  ••• 

0 

0 

0 

0 

0 

0 

0  ••• 

0 

N 

0 

0 

0 

0 

0  ••• 

-N 

0 

(30) 


where  (1  is  constructed  by  stacking  Qq  P  times  along  the  diagonal. 

Thus  we  have  2N  +  1  equations  for  the  2N  +  1  coefficients  and  the  frequency  to.  We  need 
an  additional  condition  to  solve  for  the  frequency.  Following  Bonani  and  Gilli  (1999),  we  set 
,  =  0,  which  amounts  to  saying  that  the  initial  phase  of  the  oscillator  is  arbitrary. 

We  can  solve  for  the  2N  +  1  coefficients  and  the  frequency  c o  with  the  following  set  of  at  least 
2 N  +  2  equations 


fc vQxf  —  r/(r  Lxf)  =  o 

l  xln  =  0 


(31) 


The  last  condition  amounts  to  saying  that  the  phase  of  the  oscillator  is  arbitrary  Bonani  and  Gilli 
(1999). 

For  rapid  numerical  solution,  we  can  also  write  the  Jacobian  of  the  system  Q{x.f  ,u)  from 
Eqs.  (31), 

dQ 

d  x-f  uj 


cuff  -  TDfT-1  f 2x7 

[0  0  1  0  •••  0]  0 


where  Df  is  an  M P  x  MP  matrix  with  the  following  elements  arrayed  along  the  diagonal 


9JL 

dxi 


(tl) 


8L(t 

flip  w 


M 


(33) 
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The  system  converges  pretty  well  using  Matlab’s  f  solve  function,  except  that  there  are  multi¬ 
ple  zeros  when  c o  =  (2nk)/T,  where  T  is  the  fundamental  period  and  k  >  1.  It’s  easy  to  check, 
though,  because  the  higher  harmonics  should  be  lower  in  magnitude  than  the  first  harmonics,  by 
definition,  if  u  represents  the  fundamental  period: 

xL  <  xin  and  xiik  <  xLi  (34) 

for  all  k  >  1. 


B  Harmonic  balance  technique  for  estimating  Floquet  modes 
and  multipliers 


Having  set  up  this  whole  big  framework,  it  becomes  trivial  to  estimate  both  the  Floquet  modes  and 
multipliers.  First,  we  consider  the  system  from  Eq.  (18)  and  linearize  around  the  limit  cycle,  so 
that  z (t)  =  x(t)  —  x*(f),  where  x*(f)  is  the  limit  cycle.  Then 

z  —  A(f)z  =  0,  (35) 


where  A(f)  is  the  Jacobian  of  the  state  equation  defined  on  the  limit  cycle, 


«>  -  s 


x*(t) 


If  ii/.  (f  )  is  a  Floquet  mode  of  the  system,  then 

z  (t)  =  eMfe<u  k(t) 


(36) 


(37) 


is  a  solution  of  Eq.  (35)  with  initial  condition  z(0)  =  Ufc(0).  Substituting  in  to  Eq.  (35),  we  find 
that 

fifcu k(t)  +  u k(t)  -  A(t)uk(t)  =  0.  (38) 

Then  we  again  discretize  at  M  time  intervals  to  write  u^.,  a  time-sampled  version  of  u k(t),  and 
the  corresponding  frequency  components  u{. 


u. 


—  Wfci(fi)  Ukl(t2)  ■■■  Uki(tM)  uk2  tl  ■■■  Ukp(tM) 


U 


/  _ 


=  uiiQ  uLn  uLn 


UkclN  UkslN  U20  Ukc21  Uks  21 


UkcPN  UksPN 


These  are  related  by  the  same  Y  matrix  as  before: 


u. 


=  r 


and  u{  =  Tu 


(39) 

](40) 

(41) 


We  need  to  write  a  version  of  the  A  matrix  sampled  at  M  time  points.  If  AlJ(t)  is  the  i ,  j 
component  of  A  at  time  t,  then  we  write  AT  =  diag  {Ajj(fi)  Aij(t2)  ■  ■  ■  A^-^m)},  again 

taking  each  time-sampled  value  and  aligning  them  on  the  diagonal.  The  entire  time  sampled  M P  x 
MP  matrix  Af  consists  of  the  diagonal  blocks 


'Af 

^-11 

A4 

-^-12 

A4 

^-ip 

A 4 

-^-21 

A4 

^-22 

A4 

^-2  P 

A4 

/Vpi 

A4 

AP  2 

A4 

app 

(42) 
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Having  done  that,  we  can  write  Eq.  (38)  in  the  frequency  domain 


T  [nkul  +  ujj.  -  AX] 

(43) 

/Jku[  +  uQu{  -  rA^-W* 

(44) 

4k +  cuHuf  -  Au£ 

(45) 

where  A  =  TAfT  1 .  Solving  for  the  Floquet  modes  and  multiplier  is  now  just  an  eigenvalue 
problem: 

(cuH  -  A  j  u{  =  -/ifcuf.  (46) 

There  are  (2 N  +  1  )P  eigenvalues  of  Eq.  (46),  which  are  arrayed  in  columns  along  the 
imaginary  axis  and  are  related  to  the  P  Floquet  multipliers  ///,.  of  Eq.  (35)  by  the  following 

4fk  =  /A  +  gjw  (47) 


where  q  =  ±1 . . .  N . 
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